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Abstract 

We numerically investigate the approach to the stationary state in the non- 
conservative Olami-Feder-Christensen (OFC) model for earthquakes. Starting 
from initially random configurations, we monitor the average earthquake size 
in different portions of the system as a function of time (the time is defined 
as the input energy per site in the system). We find that the process of self- 
organisation develops from the boundaries of the system and it is controlled 
by a dynamical critical exponent z ~ 1.3 that appears to be universal over a 
range of dissipation levels of the local dynamics. We show moreover that the 
transient time of the system ttr scales with system size L as ttr ~ ■ We 
argue that the (non-trivial) scaling of the transient time in the OFC model is 
associated to the establishment of long-range spatial correlations in the steady 
state. 

PACS numbers: 05.65.+b, 45.70.Ht, 89.75.Fb 

I. INTRODUCTION 

The idea of self-organised criticality (SOC) was introduced by Bak, Tang and Wiesenfeld 
(BTW) |I|] as a possible paradigm for the widespread occurrence in nature of scale free 
phenomena. It refers to the intrinsic tendency of extended, non-equilibrium systems to 
spontaneously self-organise into a dynamical critical state. In general, SOC systems are 
driven externally at a very slow rate and relax with bursts of activity (avalanches) on a 
very fast (almost instantaneous) time scale. The standard signature of SOC is a power law 
distribution of avalanche sizes and in this sense is the system said to be critical. Typical 
physical realizations of this phenomena includes, among others, earthquakes, forest fires and 
biological evolution (for a review, see e.g., Ref. 0Q]). 

A number of simple lattice models have been developed to test the applicability of SOC 
to a variety of complex interacting dynamical systems [0,0. In general these models reach a 
stationary critical state after a sufficiently long transient time. Not much attention though 
has been paid to the self-organisation process and studies have mainly concentrated on the 
properties of different models at stationarity. This is partly justified by the fact that some 
of the most studied models, such as the BTW [ffl] and the Zhang [ffl models, display a very 



1 



simple transient time behaviour. Indeed in these models the relaxation time does not scale 
with system size (time is defined as the input energy per site in the system) |^. On the 
other hand, there are also models with a more complex behaviour (see, e.g., Ref. 0). 

In this paper we investigate the approach to the stationary state in the so-called Olami- 
Feder-Christensen (OFC) model for earthquake dynamics |]^. This model has in recent 
years attracted a considerable deal of attentions especially because it has been proposed as 
an example of a system displaying SOC behaviour even with a non-conservative dynam- 
ics ||8|-[TI[|. One of the most important question in this field is indeed whether a conserving 
local dynamics is a necessary condition for SOC [^,0. For example, it is well known that 
the BTW model is subcritical if dissipation is introduced The presence of criticality in 
the non- conservative OFC model is still debated ||T5| , p!6[| . Recent numerical investigations, 
though, have shown that the model on a square lattice displays scaling behaviour, up to 
lattice sizes presently accessible by computer simulations 

The present investigation complements previous analysis of the OFC model which were 
based on the study of the probability distribution for earthquake sizes. It provides further 
numerical support in favour of criticality in the non-conservative regime. Indeed, we will 
show that the model displays a non trivial transient time behaviour: the relaxation time 
scales with system size and it is controlled by a dynamical critical exponent z that appears 
to be universal over a range of dissipation levels of the local dynamics. Moreover we will 
establish the presence of long range spatial correlations in the system. In so doing, we will be 
able to gain some insight into the mechanisms behind criticality in non-conserving systems, 
mechanisms that are very different from those at work in systems with a conservation law. 

The plan of the paper is as follows. In section II we describe the model and briefly 
summarise previous foundings relevant to our investigation. In section III we deflne the 
quantities of interest and present the results of our numerical study. Finally, in section IV 
we discuss our main conclusions. 



II. THE MODEL 

The OFC model is a coupled map lattice model, where to each site (i, j) of a square 
lattice of linear size L is associated a real variable Fij. In the initial state, at time t = 0, 
the values of the are chosen randomly in the uniform interval (0, Fc) . Subsequently the 
variables evolve according to the following two-steps dynamical rules: (i) if all sites in the 
system are stable (i.e., Fj^ < Fc), they increase simultaneously and uniformly at a constant 
rate 

(1) 

(ii) as soon as one of them reaches the threshold value Fc, the uniform driving is stopped 
and an "earthquake" starts: 

F« > ^ I ^« ^^p„ (2) 

where "nn" denotes the set of nearest neighbour sites of and a is a parameter that 
controls the level of conservation of the dynamics (a = 1/4 corresponds to the conservative 
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case). The "toppling" rule (0) can possibly produce a chain reaction, which ends when there 
are no more unstable sites in the system. At that point, the uniform growth (|l]) starts again. 
In the following we will assume, without loss of generality, a unit growth rate, i.e. v = 1. 
A crucial point in the description of the model is the choice of boundary conditions and, in 
accordance with previous investigations, we will consider open boundary conditions. These 
conditions imply that sites close to the boundaries topple according to (|) but have a smaller 
coordination number. 

There is a clear separation of time scales in the system: earthquakes occur instanta- 
neously on the slow time scale of the driving. The time in the system is therefore set by the 
slow time variable t. By construction, moreover, the time coincides with the input energy 
per site in the system. We will make use of this latter observation when we will try to 
compare the behaviour of the OFC model with other models, such as the BTW model or 
the Zhang model. 

After a sufficiently long transient time, the system settles into a stationary state, where 
the statistical properties of the model (e.g. the probability distribution for earthquake sizes) 
do not depend on time. In the BTW model and in the Zhang model, the transient time is 
relatively brief and does not scale with system size. On average, an input energy proportional 
to system size is needed to reach the steady state On the contrary, transient times in 
the OFC model are known to be extremely long, especially for large lattices. For example, 
it was claimed in Ref. |§] that 4 x 10^ earthquakes are not enough to reach stationarity in 
a system of size L = 200 for a = 0.1. This conclusion was reached by observing the very 
slow convergence of the mean earthquake size to an asymptotic value during the transient 
time. In Ref. though, a systematic investigation of the relaxation times for different a 
and L was not attempted. A more quantitative approach to the problem was proposed in 
Ref. |]TD|] by Middleton and Tang (MT). According to these authors a "self-organised" region 
develops first close to the boundaries and propagates thereafter into the bulk of the system. 
The distance from the boundaries of the invasion front grows with time as a power law, 
d{t) ~ with 7 = 0.23 ± 0.08,0.63 ± 0.08 for a = 0.07,0.15, respectively. The system 
reaches stationarity when the SOC region crosses the whole sample {d{t) ~ L). Assuming 
that the power law behaviour of d{t) holds till saturation, than the transient time of the 
system should scale as ttr ~ L^/'y("\ More recently, it has been suggested in Ref. 0] that 
two distinct relaxation times exist in the system, associated respectively with the power law 
region and the "tail" (induced by finite size effects) of the distribution of earthquake sizes. 
According to this study, the former should stabilise much faster than the latter. 



III. RESULTS 

Most of the studies on the OFC model at stationarity have concentrate on the probability 
distribution of earthquake sizes, -Pl(s), where L is the size of the system and s is the total 
number of sites that topple during an earthquake. This probability distribution does not 
show simple finite size scaling, at least in the range of lattice sizes accessible to simulations at 
present fl^. In a recent paper 0], we have focused instead on the properties of earthquakes 



confined within a fictitious subsystem of linear size A (see fig. |T]). The model is driven 
according to its usual dynamics but only those particular earthquakes that are entirely 
contained within the subsystem are counted. We have shown that if A is sufficiently smaller 
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than L the size distribution for this subset of earthquakes, Pconfi^, s), obeys ordinary finite 
size scahng, i.e. Pconf{^,s) — A~^/(s/A^), where the exponents P = 3.6 and D = 2 are 
universal over a range of values of a. 

In this work we want to address the issues briefly summarised in the previous section 
concerning the approach to stationarity in the OFC model. In order to be able to formulate 
scaling hypothesis and make use of collapse plots, we will proceed in a way similar to that 
of Ref. |TB[. We will consider earthquakes localised within given subsystems, in particular 



subsystems placed (a) at the boundaries and (b) at the centre of the system. As it is a 
prohibitive task to determine the time evolution of the entire distribution Peon/ (A, s), we 
will restrict ourselves to the mean earthquake size, < s >\^l (t), and, in general, to g-th 
moment (up to g = 4) of the distribution, < >a,l (t). To determine numerically these 
quantities we have run several simulations with different initial conditions, partitioning the 
time into bins of size At. Let n be the number of earthquakes occurring between time 
t — At/2 and t + At/2 in a given realization of the system and let Si, . . . , s„ be the sizes of 
these earthquakes. Then we define 

, < [-^1 + • • • + ^'ntf/; > . , 

< >A,L it) ^ (3) 

where < . . . > denotes an average over different realizations of the system, i.e. over different 
initial conditions. For each system, the parameter At has been choose small with respect to 
the transient time ttr but large enough to collect reasonably accurate statistics. 

We consider first the case of subsystems placed adjacent to a boundary of the system, 
in a symmetric position with respect to the corners. In figure 2 we report < s >a,l as a 
function of time for a = 0.18 and some A and L. We observe that if the linear dimension A 
of the box is sufficiently smaller than the linear dimension of the system L (approximately 
L > 4A) then the curve < s >x,l (t) becomes indistinguishable for different L. This has 
been verified also for other values of a and A. We will therefore denote with < s >a {t) the 
mean earthquake size in this limit. It is already visible from figure 2 that the relaxation time 
of < s >A (t) increases with A. This is an indication in support of the scenario proposed by 
MT. Regions close to the boundaries reach stationarity sooner, signalling that an invasion 
front is moving toward the bulk of the system. Some of MT conclusions nonetheless will 
have to be modified as we will show later. 

In order to describe quantitatively the invasion from the boundaries of the self-organised 
region we make the following simple scaling hypothesis 

< s >A (t) = X'^Fit/y) (4) 

where rj and z are two suitable critical exponents. In particular 2; is a dynamical critical 
exponent that should satisfy z = l/'-f, where 7 is the "invasion" exponent as defined by MT. 
In the limit of t ^ 00 the scaling function F{x) saturates to a constant, implying that the 
exponent rj is related to the finite size exponents of the probability distribution Pconf{X, s) 
by the relation rj = 2D — (3 0.4. In figure |^ we report collapse plots of the form (D for 
various values of a. We observe that a reasonably good collapse could be obtained for all 
the a if we choose the universal exponent z = 1.3 ± 0.1. We are therefore led to conclude 
that a universal exponent z exist contrary to the claims by MT. 
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We consider next subsystems of different sizes placed at the centre of tlie system. We 
observe that the relaxation time does not depend on the size of the subsystems (see fig. ^). 
This confirm that the self-organisation mechanism develops from the boundaries and that the 
system enters stationarity when the self-organised region span the whole system. Only when 
the bulk of the system is reached by the self-organised region is stationarity settled so that 
concentric subsystem of different sizes will inevitably reach stationarity at the same time. 
We have also tested the relaxation times for higher moments of the avalanche probability 
distribution to see whether different parts of the distribution (e.g. power law part and 
the "tail") have different relaxation times as suggested in Ref. [|TT|. In our investigation 
though we have not observed any significant difference in the relaxation times associated 
with different moments. We report as an example in fig. the behaviour of the first, second 
and fourth moment in a particular case. 

The scaling equation (^) suggests that the transient time in the OFC model scales with 
system size as tt^ ~ L^- One way to test this is by comparing the time behaviour of the 
average earthquake size in a central subsystem of size A for different system sizes L. Indeed 
the asymptotic value < s >a,l (t oo) should not depend on L. We report as example 
the case for a = 0.18 in fig. where we have rescaled the time by a factor L^. The curves 
show some noisy behaviour, due to the difficulties in collecting good statistics (relatively few 
earthquakes occur in the bulk of the system as compared to the boundaries). Nonetheless 
the value deduced for the exponent, z ^ 1.3, is consistent with the determination made 
through the analysis of the earthquakes occurring at the boundaries. We have obtained 
similar results also for other values of a. In addition, besides the average earthquake size, 
we have considered also the time behaviour of other quantities such as the roughness of the 
energy landscape (in analogy to surface growth problems) and the number of earthquakes 
per unit time. All these different quantities on average reach stationary at the same time. 

The algebraic divergence of the relaxation time with system size reflects the presence 
of long-range spatial correlations in the stationary state. Indeed if correlations were only 
short range, than one would expect that the transient time would not scale with system 
size. This is for example the case for the BTW model in d dimensions, where the height- 
height correlations are algebraic but decays as fast as r"^*^ (r being the distance between 
two sites) fl^. We have measured for various L and a the probability distribution of the 
spatially averaged force in the system 

In a system with sufficiently short range correlations, this probability distribution would 
tend, in the limit L — > oo, to a gaussian distribution around the mean due to the central 
limit theorem (this is indeed what results in the BTW model). Let < M > and ctm be 
respectively the average and the standard deviation of the distribution. In figure |^ we have 
plotted log((jMP(M)) versus (M— < M >)/aM for various L and a. Using these coordinates 
a gaussian function would result in an inverse parabola. For each a value the data of figure 
^ collapse on a single function, which is clearly not gaussian (deviations from gaussianity are 
more pronounced for increasing a values). This indicates that the central limit theorem does 
not hold in this case, not even for large L, suggesting that long range algebraic correlations 
are present and therefore the sum (^ can not be decomposed into a sum of independent 
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terms. This observation is in agreement with the results reported in Ref. |T3[ where the 
presence of long-range spatial correlations were deduced from the behaviour of a suitably 
defined susceptibility, x = {LaMY- 1^ claimed that x diverges as and correspondingly 
that (Tm is, to leading order, independent on L (if M was a sum of uncorrelated variables, 
o"A/ would decrease as In our investigation we have found that au slightly increases 

with L but asymptotically tends to a constant value, in accordance with Ref. [0 (M is a 
bounded variable so (Jm can not grow indefinitely). 



IV. CONCLUSIONS 

In conclusion, in this paper we have examined the process of self-organisation in the 
OFC model. By considering earthquakes confined within a given subsystem we have been 
able to clarify some of the issues related to this problem. In accordance with Middleton and 
Tang |]TD[ we have found that SOC develops first close to the boundaries and subsequently 
invades the interior of the system. The invasion process is controlled by a dynamical critical 
exponent, z ~ 1.3, which, contrary to previous claims, is universal over a range of values of 
the dissipation level of the local dynamics. We have shown moreover that the transient time 
in the system scales with system size as ttr ~ ■ This is a peculiarity of the OFC model 
as other "sandpile-like" models (e.g. the BTW and the Zhang models) do not display any 
scaling in the transient time. We have associated this feature with the presence of long-range 
spatial correlations in the stationary state. 



Our findings are in general agreement with recent works on the OFC model |17| , p!8 
Indeed we have provided complementary evidences (not based on the probability distribution 
for earthquake sizes) that the model is critical even in a non conservative regime. Moreover it 
confirms that there is universality in the system and that finite-size scaling can be recovered 
by considering subsystems whose linear extent is sufficiently small. 

Finally, it is interesting to remark that the probability distribution for the spatially 
averaged force in the system is somewhat reminiscent of a probability distribution observed 
in a confined turbulent fiow experiment (HHP). As a term of comparison we have 
reported in fig. ^the BHP functional form over-imposed to the curve for a = 0.21. Attempts 
to link SOC systems to turbulent phenomena have long been suggested, but only recently 
this has been put on a firmer basis ||21||. 
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FIGURES 
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FIG. 1. Schematic representation of earthquakes entirely confined within a subsystem of hnear 
size A (dashed Une). ToppUng sites are denoted with a cross. 
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FIG. 2. Average earthquake size < s >a,l in a subsystem placed at the boundary as a function 
of time; a = 0.18 and, from bottom to top, A = 32, A = 64 and A = 128. 
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FIG. 3. Collapse plots of < s >\ (t) for a subsystem placed at the boundary and for (a) 
= 0.15, (b) a = 0.18 and (c) a = 0.21. The value of the dynamical critical exponent is z = 1.3. 
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FIG. 4. Time dependence of the average earthquake sizes in subsystems placed at the centre 
of the system for a = 0.18 and L = 256; (a) average earthquake size < s >a,l for, from bottom to 
top, A = 16, 32, 64, 128, 256 and (b) q — th moment of the distribution for A = 32 and, from bottom 
to top, g = 1, 2,4. 
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FIG. 5. Collapse plots of < s >a,l (*) for a subsystem of size A = 32 placed at the centre of 
a system of size L; the conservation parameter is a = 0.18. The value of the dynamical critical 
exponent is z = 1.3. 
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FIG. 6. Rescaled probability distribution of the spatially averaged force in the system for, 
from bottom to top, a = 0.21,018,0.15. < M > and fTj\/ are respectively the average and the 
standard deviation of the distribution. The top and bottom curves have been shifted by a factor 
of 10 respectively up and down for visual clarity. Squares represent the BHP (rescaled) probability 
distribution observed in an experiment on turbolence [^0| . 
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